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We design an irreversible Markov chain Monte Carlo algorithm for the self-avoiding walk (SAW). It outper- 
forms the Berretti-Sokal algorithm. The gained efficiency increases with the spatial dimension, from about 
10 times in two dimensions to around 40 times in five dimensions. The algorithm violates the widely used 
detailed balance condition and satisfies the weaker balance condition. We employ the irreversible method to 
study the finite-size scaling of SAW above the upper critical dimension. 
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|. INTRODUCTION 


The self-avoiding walk (SAW) serves as an important 
model in polymer physics (see e.g. Ref. [il] and references 
therein). In the grand-canonical ensemble, the walks can 
have different lengths and the model is defined by the 


partition sum 
c= Sa", (1) 
w 


where the summation is over all possible self-avoiding 
paths, with |w| being the length of a walk w. It is equiv- 
alent to the n + 0 limit of the O(n) mode. Markov 
chain Monte Carlo (MCMC) methods have been exten- 
sively used in the simulation of SAW. 

There are two key ingredients for designing a MCMC 
algorithm. One is the balance condition (BC) which tells 
that the probability flow entering into a configuration 
equals the flow out of the configuration, thus it is guaran- 
teed that the chain converges to a stationary distribution. 
And the other one is the ergodicity which requires that all 
configuration states can be reached from a random initial 
configuration. In practice, the BC is usually satisfied by 
employing the detailed balance condition (DBC), which 
means that the probability flow from a configuration to 
another one equals the reverse flow, i.e. the dynamics is 
reversible. The DBC is sufficient but not necessary. 

Recently there appear several successful studiest+4 
which show promising future for MCMC algorithms be- 
yond the DBC: geometric allocation approaches have 
been applied to the Potts model4243; irreversible 
MCMC methods have been designed for the mean- 
field Ising mode; event-chain Monte Carlo (ECMC) 
methods have been proposed for simulation of hard- 
sphere systems®£ and generalized to particle systems 
with arbitrary pairwise interactions®, including the soft- 
disk systems2, the XY modclH2 and the Heisenberg 
mode. The geometric allocation method outperforms 
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the Metropolis-Hasting (MH) method by 6.4 times for 
the q = 4 square lattice Potts model, and is increasingly 
advantaged as q increases#2. For the mean-field Ising 
model, the irreversible MCMC methods have a dynamic 
exponent z ~ 0.85, which is much smaller than z ~ 1.43 
for the reversible MH method42. Comparing with the 
MH method, the speedup of the ECMC method reaches 
two orders of magnitude in large systems consisting of 
10° hard spheres£. And for the three-dimensional fer- 
romagnetic Heisenberg model, it has been reported that 
the ECMC method has a dynamic exponent z ~ 1, in 
contrast to z ~ 2 for the MH methodt. 


While local probability flow circles are introduced in 
the geometric allocation approach, much larger or even 
global circles appears in other methods mentioned above. 
This is achieved by making use of the lifting technique. 
This technique makes replications of the system. In each 
replica, the probability flow goes either in a chosen direc- 
tion or to another replica. There exists net probability 
flow in each replica, and between different replicas, but 
the combined flow is balanced, i.e. the BC is satisfied. 
For example, in simulating the mean-field Ising model4, 
two replicas are used. In one replica, only the positive 
spins are proposed to flip (the total magnetization M al- 
ways decreases), and only the negative spins are proposed 
to flip in the other replica (M always increases). If the 
spin-flip update is rejected, switching of replica will be 
proposed. Updates in one replica can persist very far un- 
til a spin-flip proposal is rejected, after which the replica 
changes and updates can persist in the other replica un- 
til another spin-flip is rejected. Thus large circles are 
introduced in the probability flow. In this way, the diffu- 
sive feature in phase space is suppressed or even replaced 
by ballistic-like behavior. And due to this fact, for the 
mean-field Ising model, comparing with the MH method, 
the dynamics of the lifted Markov chain method is qual- 
itatively faster, i.e. the dynamic exponent is z ~ 0.85 
instead of z ~ 1.43 for the MH method42. However, it 
is found that the lifting trick does not help much for the 


Ising model in low dimensions244. 


In this work, we design an efficient MCMC algorithm 
without DBC for SAW by employing the lifting tech- 


nique with two replicas. The change of the SAW path 
is different in the two replicas: it cannot increase its 
length in one replica and cannot decrease its length in 
the other. And switching of replicas is allowed to enable 
formation of probability flow circles. This new algorithm 
is irreversible since it violates the DBC. Numerical results 
show that the irreversible MCMC method is far superior 
to the Berretti-Sokal (BS) algorithm+£. We use this new 
method to explore the finite-size scaling (FSS) of SAW 
above the upper critical dimension. 

The remaining part of this article is organized as fol- 
lows. In Sec. [L| we recall traditional Monte Carlo meth- 
ods for SAW, including the MH algorithm and the BS 
algorithm. These algorithms preserve reversibility as in- 
dicated by the DBC. Section [M] describes in detail our 
irreversible MCMC algorithm for SAW. We compare the 
performance of the above algorithms in Sec. [IV] Section[V] 
contains a FSS analysis for SAW on a periodic hypercubic 
lattice in five dimensions, which is above the upper crit- 
ical dimension de = 4. A brief conclusion and discussion 
is presented in Sec. 


Il. REVERSIBLE MCMC ALGORITHMS 


We introduce below the MH and BS algorithms for 
SAW. These two algorithms employ the DBC, thus their 
dynamics are reversible. 


A. Metropolis-Hasting algorithm 


Let # be a regular d-dimensional lattice with coordi- 
nation number z. Then an N-step SAW w on Z (N = 
|w|) is a sequence of distinct sites w(0), w(1), ... ,w(N) 
in Y such that each site is a nearest neighbor of its pre- 
decessor. We assume that all walks have an open end 
fixed at the origin, i.e. w(0) = 0, and the other open 
end denoted by J. The MH algorithm+248 is a standard 
reversible MCMC algorithm, which is constructed as be- 
low. 


MH Algorithm. 

(1) Randomly choose one of the z neighboring sites of J, 
say I’, and propose a symmetric update to for the edge 
connecting J and J’, i.e. propose moving the open end I 
to I’. 


(2) Accept the proposal with probability 
min{1l,7(w’)/7(w)}, where w’, w represent the 
two walks with open end I’ and I, respectively, 


and m(w) gives the weight of a configuration w. 


Here the weight 7 is given by 


all 
rw) = xu), (2) 


where 


ne 1 ifw isan SAW (3) 
XW) =) 0 if w is not an SAW. 


The acceptance probabilities are given by 

P(AN = +1) = min{1, r}y(w’) (4) 
and 

P(AN = —1) = min{1,1/z}, (5) 


where AN = |w’| — |w|. The a priori probability A(w > 
w’) describes the probability of proposing a transition 
from configuration w to w’. Here we have A(w > w’) = 
1/z for both AN = +1 and AN = —1. It can be verified 
that the DBC is satisfied 


pw > w’) = ou’ > w) (6) 


where ¢(w > w) = m(w)A(w > w’)P(w —> w’) is the 
local probability flow. 


B. Berretti-Sokal algorithm 


An efficient MCMC algorithm for the SAW is the BS 
method42. The BS algorithm uses a priori probabilities 
different from the MH method. 


BS Algorithm. 

(1) Propose randomly to increase the walk length 
(AN = +1) or to decrease it (AN = —1), with equal 
probability. 

(2a) When AN = +1 is proposed, one randomly selects 
a neighboring site of J, say I’. If the proposed move 
leads to self-intersection, it is rejected and the old 
configuration is counted again in sampling (a “null tran- 
sition”). Otherwise, with probability Pgs(AN = +1), 
accept the proposal and move the open end from I to I’. 
(2b) When AN = —1 is proposed, if the proposal is 
made to delete a bond from an already-empty walk, it 
is rejected and a “null transition” presents. Otherwise, 
with probability Pgs(AN = —1), accept the proposal, 
delete the last bond and move the open end from I to 
its predecessor I’. 


The predecessor of site I is excluded from choice in (2a), 
thus each I’ is selected with probability 1/(z — 1). The 
a priori probabilities are A(w —> w’) = 1/2(z — 1) for 
AN = +1, and A(w > w’) = 1/2 for AN = —1. The 
acceptance probabilities are chosen as 


Pes(AN = +1) = min{1, z(z — 1)} (7) 
and 
Pgs(AN = —1) = min{1, 1/z(z — 1)}. (8) 


It can be verified that the above A and Pgs also satisfy 
the DBC given in Eq. (6). 


We note that when M = 0, there are z candidates for 
I’, instead of z — 1. A convenient treatment is to set for 
N = 0 one of the z neighbors not available. It is also 
noted that there is a slight difference between the BS al- 
gorithm described above and the original BS algorithm. 
In the original algorithm, a priori probabilities are cho- 
sen such that acceptance probabilities always equal one 
if the proposed step does not lead to self intersection 
or does not delete a bond from an already-empty walk. 
However, since we focus on properties near the critical 
point xe, which has its value very close to 1/(z — 1), this 
does not lead to significant difference in the efficiency. 


Ill. IRREVERSIBLE MCMC METHOD 


The BS algorithm in the previous section clearly de- 
compose the proposed steps into those with AN = +1 
and others with AN = —1. To get an irreversible algo- 
rithm, a direct idea is to use two replicas, labeled with 
(+) and (—), and only propose to increase the walk length 
in the (+) replica, and to decrease the walk length in the 
(—) replica. The DBC is violated since in each replica 
the probability flow only goes in one direction. To en- 
sure convergence to the equilibrium distribution, we use 
the BC instead, which is achieved by allowing transitions 
between the two replicas. Figure[I]shows the probability 
flow. The irreversible algorithm goes as follows. 


Irreversible Algorithm. 

(1) Propose to increase the walk length AN = +1. 
Randomly select a neighboring site of I, say I’. If the 
proposed move leads to self-intersection, it is rejected. 
Otherwise, with probability Pt, accept the proposal 
and move the open end from J to I’. Once the proposal 
is rejected, the walk switches to the (—) replica. 

(2) Propose to decrease the walk length AN = —1. If 
the proposal is made to delete a bond from an already- 
empty walk, it is rejected. Otherwise, with probability 
P~, accept the proposal, delete the last bond and move 
the open end from IJ to its predecessor I’. Once the 
proposal is rejected, the walk switches to the (+) replica. 


To observe the probability flow, we count the flow into 
and out of a configuration in the (+) replica with length 
N. The input flow from a M — 1 configuration is 

1 
Z= 


ot (NM -1 >N) =^! PH. (9) 


There also exists an input flow from the M-step walk in 
the (—) replica 
EINT >N?) = (1 - P7). (10) 


The outgoing flow decomposes into two parts, one part 
is that leads to SAWs with length M + 1 


Ptal N >N +1) = rN —— 


t 
ou T 


Pr, (11) 


FIG. 1. Probability flow of the irreversible Monte Carlo algo- 
rithm for SAW. 


where 2’ (0 < z’ < z — 1) is the number of valid (M + 1)- 
step SAWs that can be obtained from the present M-step 
SAW. The other part is that leads to the corresponding 
N-step walk in the (—) replica, due to rejection of the 
proposed AN’ = +1 walk in the (+) replica, 


dour (NF >N) 


/ = if 
Z (1 Bae = 
z—1 


where MT and M+ denote the N-step walk in the (—) 
and (+) replica, respectively. 
The BC condition implies that 
PEN -13N)+¢5.WN7 >N?) 


= Gla N >N 1) + Plua NT > NT). (13) 


The above equation is satisfied if we choose Pt = 
Pes(AN = +1) = min{1,z(z — 1)} and P% = 
Pgs(AN = —1) = min{1,1/æ(z — 1)}. Similarly, one 
can verify that the BC is satisfied for a M-step walk in 
the (—) replica. The above irreversible algorithm is easy 
to implement — only small changes to the BS code are 
needed, due to the similarity in their descriptions. 

We note that this irreversible algorithm satisfies the 
skew DBG48, which requires that the probability flow 
from a configuration to another in one replica equals the 
reverse flow in the other replica, i.e. 6T(NM > N +1) = 
eg (N+13N). 


=^ 


(12) 


IV. PERFORMANCE 


We compare the efficiency of the algorithms by observ- 
ing the integrated autocorrelation time 7, which is defined 


byl 
esc serie (2 -0’). (14) 
2 


Here ôO represents the statistical error for an arbitrary 
observable O, At is the time interval at which the sam- 
ples are taken, and n = tmazx/At is the total number 
of samples, where tmar is maximum measure time. The 
quantity 7 indicates how many Monte Carlo steps are 
needed to generate an effectively independent sample. 


TABLE I. Value of critical point x, in different spatial dimen- 
sions d. 


d Xe 
2 0.379 052 277 758(4)22 
3 0.213 491 0(3)24 
4 0.147 622 3(1)}3% 
5 0.113 140 81(4)22 
0.113 140 84(1) [this work] 
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FIG. 2. Autocorrelation time of the Metropolis-Hasting (MH) 
algorithm rmn divided by that of the Berretti-Sokal (BS) or 
irreversible (IR) algorithm, versus the linear system size L. 
The autocorrelation time is calculated from observing the av- 
erage walk length N = (NV), at the critical point £e, for two 
to five dimensions. For each linear size L, we ran 20 indepen- 
dent jobs. In each job, 50 x 1024 x L“/2% steps were thrown 
for thermalization, and another 200 x 1024 x L? / 2¢ steps were 
accounted for sampling. 


We conducted simulations at the critical point ze (see 
Table [I] for their values; simulations for d = 5 were con- 
ducted at £e = 0.113 140 85.) for the square, cubic and 
d-dimensional hypercubic lattice with the spatial dimen- 
sion d = 4,5, using periodic boundary conditions. 

Figure 2] compares the autocorrelation time 7 for dif- 
ferent algorithms, calculated from observing the average 
walk length N = (NV). It can be seen that the perfor- 
mance of the irreversible algorithm is far superior to the 
BS and the MH algorithms. The efficiency of the irre- 
versible algorithm increases as the spatial dimension d 
becomes larger. For d = 2, it outperforms the MH al- 
gorithms about 12 times, and outperforms the BS algo- 
rithm about 8 times. And for d = 5, it outperforms the 
MH algorithms over 200 times, and outperforms the BS 
algorithm around 40 times. 

We also compare the efficiency of different algorithms 
in units of CPU time using the same machine (a desktop 
computer with 3.8 GiB memory and four Intel i5 cores). 
The CPU time consumed in T sweeps with £x = £e are 
shown in Table [II] for different algorithms. The results 


TABLE II. CPU time (in units of seconds) needed in 7 sweeps 
for different algorithms at x = ze, measured in d-dimensional 
systems with linear size L. For convenience, the ratio of CPU 
time tps/tir is calculated, as well as the ratio of autocorrela- 
tion time Tgs/TIR- 


dL tr tes tMH tgs/tiR TBs/TIR 
2 1024 0.130(11) 2.0(3)  3.3(3) 15(2) 8(2) 
3 128 0.014 8(7) 0.264(11) 0.81(14) 18(2) 13.2(11) 
464 0.043(2) 0.73(3) 2.5(3) 17(2)  24(2) 
5 32 0.035(2) 1.21(12) 4.3(9) 35(6) 43(6) 
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FIG. 3. Autocorrelation time of the Metropolis-Hasting (MH) 
algorithm 7 divided by that of the Berretti-Sokal (BS) or 
irreversible (IR) algorithm, versus the linear systems size L. 
The autocorrelation time is calculated from observing the re- 
turning probability Do = (Do), at the critical point xe, for 
two to five dimensions. 


are similar to those from the comparison of the autocor- 
relation time. 

From Eq. (14), different observables may lead to differ- 
ent estimations of the autocorrelation time 7. The results 
above are obtained from observing the walk length. In 
order to account for this difference, we measure another 
quantity, i.e. the returning probability Do = (Do) which 
is defined as the probability that the walk returns to its 
start point (thus having zero walk length). The results 
are shown in Fig.[3] It can be seen that the improvement 
of efficiency is similar to those obtained from observing 
the walk length. 


V. FINITE-SIZE SCALING ABOVE THE UPPER 
CRITICAL DIMENSION 


While the FSS of SAW is usually presented in terms of 
the walk length or the end-to-end distance, here we an- 
alyze the dependence of observables on the linear size of 
the lattice L. An example on this kind of FSS has been 
presented for three-dimensional SAW in Ref. p3, where it 
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FIG. 4.  Finite-size scaling near the critical point £e = 
0.113 140 84. The left panel shows the rescaled average walk 
length N/L°/? versus x (top), and the specific heat C versus 
x (bottom). The right panel shows the same quantities using 
a rescaled horizontal coordinate (x — £e) L5, data points at 
different (x, L) collapse well into a single scaling function as 
L becomes large. 


is shown that dimensionless ratios nicely intersect at the 
critical point ze. Here we use the irreversible method 
to investigate the FSS above the upper critical dimen- 
sion de = 4. For the Ising model with periodic boundary 
conditions, Binder et al24 predicts that above the up- 
per critical dimension, the thermal and magnetic scaling 
exponent are y = d/2 and yp = 3d/4 respectively. We 
demonstrate below that the same exponents also apply 
to SAW. 

Using the irreversible method, extensive simulations 
were done for SAW in five dimensional hypercubic lattice 
with periodic boundary conditions. The average walk 
length N is an energy-like quantity, which scales as 


2 
N = L” (ao +X ape" L +b L”), (15) 
k=1 
where € = g% — ze and yı is the leading correction ex- 
ponent, Go, @1, @2, bı are non-universal constants. We 
also sampled the specific-heat C = L~4 ((N?) — (N)?), 
which behaves as following 


2 
C= LY 4 (ag + So apt L +b L”). (16) 
k=1 

The MC data were fitted to Eqs.(I5) and (16), with dif- 
ferent values of yı, leading to ze = 0.113 140 84(1) and 
values of y, consistent with 5/2. The value of ze agrees 
with the best previous estimates £e = 0.113 140 81(4)22. 
Figure [4] shows that the data collapse well onto a single 

scaling function when plotted in rescaled coordinates. 
We also measured the returning probability Do at the 


critical point, which is expected to scale as22 


Do = L2 2Yh (ao + b L” + ate .) ; (17) 
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FIG. 5. Plot of DoL?¥~% vs L at the critical point te = 
0.113 140 85, with yn = 3d/4, d = 5. 


where yp = 3d/4 = 15/4, and higher-order correction 
terms are omitted. The above scaling behavior is demon- 
strated in Fig. B] where DoL?¥"~¢ is plotted versus L. 
The fact that DoL?¥"—4¢ converges to a nonzero value as 
system size increases indicates that the leading scaling 
behavior of Do is indeed ~ L¢~2¥%, 

The irreversible method can also be used to simulate 
SAW on lattices with free boundary conditions. Above 
the upper critical dimension, the above values y, = d/2 
and yn = 3d/4 account for effects of the dangerous irrel- 
evant field24, otherwise the exponents should be y} = 2 
and y, = d/2+1. The latter values have been observed 
for SAW on five-dimensional hypercubic lattices with free 


boundary conditions2%. 


VI. CONCLUSION AND DISCUSSION 


We construct an irreversible Monte Carlo algorithm 
for SAW. It violates the DBC, and satisfies the weaker 
BC. Comparing with the MH and BS algorithms, both 
being reversible, the irreversible method is much more ef- 
ficient. While the BS algorithm is about d times more ef- 
ficient than the MH algorithm, the irreversible algorithm 
is much more superior to the BS algorithm. Comparing 
the irreversible method with the BS method, the gain 
of efficiency increases with the spatial dimension: from 
around 10 times in 2D to around 40 times in 5D. 

The irreversible algorithm introduces an auxiliary vari- 
able to enlarge the configuration space, and as a re- 
turn, this gives one some freedom in designing MCMC. 
More specifically, the irreversible algorithm makes use 
of two replicas. While the BS algorithm proposes ei- 
ther the AN = +1 move or the AN = —1 move with 
a nonzero probability, the irreversible algorithm pro- 
poses only AN’ = +1 move in one replica, and only the 
AN = —1 move in the other replica. The probability 
flow in each replica is irreversible. In order to satisfy the 
BC, switching of replica is employed when a proposed 
AN = +1 or — 1 move is rejected. The implementation 
of the irreversible algorithm only needs small modifica- 


tions to the BS code. It is expected that the irreversible 
MCMC method may be useful in studies of interacting 
SAW, other polymeric systems and soft matter. 
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